!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine tvd_lf (ro, mx, my, mz, bx, by, bz, en, nx, ng, flux_bnd, n_neighs, sys, dx, dt, gm, riemann_solver_flag)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ng, n_neighs, sys, riemann_solver_flag
   real(8), intent(in) :: dx, dt, gm
   real(8), dimension(n_neighs, sys), intent(inout) :: flux_bnd
   real(8), dimension(nx), intent(inout) :: ro, mx, my, mz, bx, by, bz, en

   integer(4) :: i, ip1, im1

   real(8), dimension(nx) :: ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre
   real(8), dimension(nx) :: ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l
   real(8), dimension(nx) :: ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r
   real(8), dimension(nx) :: fro, fmx, fmy, fmz, fbx, fby, fbz, fen

   real(8) :: fro_l, fmx_l, fmy_l, fmz_l, fbx_l, fby_l, fbz_l, fen_l
   real(8) :: fro_r, fmx_r, fmy_r, fmz_r, fbx_r, fby_r, fbz_r, fen_r
   real(8) :: limiter
   integer(4) :: sgn
   real(8) :: dtdx, u1, u2, vx, vy, vz, b2, v2, c2, s2, ca2, cfx, pr, max_speed, gmm1

!----------------------------------------------------------------------------------------------------------------------------------|
   dtdx = dt / dx
   gmm1 = gm - 1.0d0

   do i = 2, nx - 1
      ip1 = i + 1
      im1 = i - 1

      u1 = ro(i) - ro(im1)
      u2 = ro(ip1) - ro(i)
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (abs (u1), sgn * u2))
      ro_l(i) = ro(i) + 0.5d0 * limiter
      ro_r(i) = ro(i) - 0.5d0 * limiter

      u1 = mx(i) - mx(im1)
      u2 = mx(ip1) - mx(i)
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (abs (u1), sgn * u2))
      mx_l(i) = mx(i) + 0.5d0 * limiter
      mx_r(i) = mx(i) - 0.5d0 * limiter

      u1 = my(i) - my(im1)
      u2 = my(ip1) - my(i)
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (abs (u1), sgn * u2))
      my_l(i) = my(i) + 0.5d0 * limiter
      my_r(i) = my(i) - 0.5d0 * limiter

      u1 = mz(i) - mz(im1)
      u2 = mz(ip1) - mz(i)
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (abs (u1), sgn * u2))
      mz_l(i) = mz(i) + 0.5d0 * limiter
      mz_r(i) = mz(i) - 0.5d0 * limiter

      u1 = bx(i) - bx(im1)
      u2 = bx(ip1) - bx(i)
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (abs (u1), sgn * u2))
      bx_l(i) = bx(i) + 0.5d0 * limiter
      bx_r(i) = bx(i) - 0.5d0 * limiter

      u1 = by(i) - by(im1)
      u2 = by(ip1) - by(i)
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (abs (u1), sgn * u2))
      by_l(i) = by(i) + 0.5d0 * limiter
      by_r(i) = by(i) - 0.5d0 * limiter

      u1 = bz(i) - bz(im1)
      u2 = bz(ip1) - bz(i)
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (abs (u1), sgn * u2))
      bz_l(i) = bz(i) + 0.5d0 * limiter
      bz_r(i) = bz(i) - 0.5d0 * limiter

      u1 = en(i) - en(im1)
      u2 = en(ip1) - en(i)
      if (u1 .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif
      limiter = sgn * max (0.0d0, min (abs (u1), sgn * u2))
      en_l(i) = en(i) + 0.5d0 * limiter
      en_r(i) = en(i) - 0.5d0 * limiter
   enddo

   do i = 1, nx
      vx = mx_l(i) / ro_l(i)
      vy = my_l(i) / ro_l(i)
      vz = mz_l(i) / ro_l(i)
      b2 = bx_l(i) * bx_l(i) + by_l(i) * by_l(i) + bz_l(i) * bz_l(i)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_l(i) - v2 * ro_l(i) / 2.0d0 - b2 / 2.0d0) * gmm1
      fro_l = mx_l(i)
      fmx_l = vx * mx_l(i) + pr + b2 / 2.0d0 - bx_l(i) * bx_l(i)
      fmy_l = vx * my_l(i) - bx_l(i) * by_l(i)
      fmz_l = vx * mz_l(i) - bx_l(i) * bz_l(i)
      fbx_l = 0.0d0
      fby_l = vx * by_l(i) - bx_l(i) * vy
      fbz_l = vx * bz_l(i) - bx_l(i) * vz
      fen_l = (en_l(i) + pr + b2 / 2.0d0) * vx - (bx_l(i) * vx + by_l(i) * vy + bz_l(i) * vz) * bx_l(i)

      vx = mx_r(i) / ro_r(i)
      vy = my_r(i) / ro_r(i)
      vz = mz_r(i) / ro_r(i)
      b2 = bx_r(i) * bx_r(i) + by_r(i) * by_r(i) + bz_r(i) * bz_r(i)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_r(i) - v2 * ro_r(i) / 2.0d0 - b2 / 2.0d0) * gmm1
      fro_r = mx_r(i)
      fmx_r = vx * mx_r(i) + pr + b2 / 2.0d0 - bx_r(i) * bx_r(i)
      fmy_r = vx * my_r(i) - bx_r(i) * by_r(i)
      fmz_r = vx * mz_r(i) - bx_r(i) * bz_r(i)
      fbx_r = 0.0d0
      fby_r = vx * by_r(i) - bx_r(i) * vy
      fbz_r = vx * bz_r(i) - bx_r(i) * vz
      fen_r = (en_r(i) + pr + b2 / 2.0d0) * vx - (bx_r(i) * vx + by_r(i) * vy + bz_r(i) * vz) * bx_r(i)

      ro_pre(i) = ro(i) - 0.5d0 * dtdx * (fro_l - fro_r)
      mx_pre(i) = mx(i) - 0.5d0 * dtdx * (fmx_l - fmx_r)
      my_pre(i) = my(i) - 0.5d0 * dtdx * (fmy_l - fmy_r)
      mz_pre(i) = mz(i) - 0.5d0 * dtdx * (fmz_l - fmz_r)
      bx_pre(i) = bx(i) - 0.5d0 * dtdx * (fbx_l - fbx_r)
      by_pre(i) = by(i) - 0.5d0 * dtdx * (fby_l - fby_r)
      bz_pre(i) = bz(i) - 0.5d0 * dtdx * (fbz_l - fbz_r)
      en_pre(i) = en(i) - 0.5d0 * dtdx * (fen_l - fen_r)
   enddo

   do i = 1, nx
      ro_l(i) = ro_pre(i) + ro_l(i) - ro(i)
      ro_r(i) = ro_pre(i) + ro_r(i) - ro(i)
      mx_l(i) = mx_pre(i) + mx_l(i) - mx(i)
      mx_r(i) = mx_pre(i) + mx_r(i) - mx(i)
      my_l(i) = my_pre(i) + my_l(i) - my(i)
      my_r(i) = my_pre(i) + my_r(i) - my(i)
      mz_l(i) = mz_pre(i) + mz_l(i) - mz(i)
      mz_r(i) = mz_pre(i) + mz_r(i) - mz(i)
      bx_l(i) = bx_pre(i) + bx_l(i) - bx(i)
      bx_r(i) = bx_pre(i) + bx_r(i) - bx(i)
      by_l(i) = by_pre(i) + by_l(i) - by(i)
      by_r(i) = by_pre(i) + by_r(i) - by(i)
      bz_l(i) = bz_pre(i) + bz_l(i) - bz(i)
      bz_r(i) = bz_pre(i) + bz_r(i) - bz(i)
      en_l(i) = en_pre(i) + en_l(i) - en(i)
      en_r(i) = en_pre(i) + en_r(i) - en(i)
   enddo

   if (riemann_solver_flag .eq. 1) then
      call hlld (fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                           &
      ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, nx, gm)
   elseif (riemann_solver_flag .eq. 2) then
      call hllc (fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                           &
      ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, nx, gm)
   elseif (riemann_solver_flag .eq. 3) then
      call roe (fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                            &
      ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, nx, gm)
   else
      do i = 1, nx - 1
         ip1 = i + 1

         vx = mx_l(i) / ro_l(i)
         vy = my_l(i) / ro_l(i)
         vz = mz_l(i) / ro_l(i)
         b2 = bx_l(i) * bx_l(i) + by_l(i) * by_l(i) + bz_l(i) * bz_l(i)
         v2 = vx * vx + vy * vy + vz * vz
         pr = (en_l(i) - v2 * ro_l(i) / 2.0d0 - b2 / 2.0d0) * gmm1
         c2 = gm * pr
         s2 = c2 + b2
         ca2 = bx_l(i) * bx_l(i)
         cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_l(i) / 2.0d0)
         max_speed = cfx + abs (vx)

         fro_l = mx_l(i)
         fmx_l = vx * mx_l(i) + pr + b2 / 2.0d0 - bx_l(i) * bx_l(i)
         fmy_l = vx * my_l(i) - bx_l(i) * by_l(i)
         fmz_l = vx * mz_l(i) - bx_l(i) * bz_l(i)
         fbx_l = 0.0d0
         fby_l = vx * by_l(i) - bx_l(i) * vy
         fbz_l = vx * bz_l(i) - bx_l(i) * vz
         fen_l = (en_l(i) + pr + b2 / 2.0d0) * vx - (bx_l(i) * vx + by_l(i) * vy + bz_l(i) * vz) * bx_l(i)

         vx = mx_r(ip1) / ro_r(ip1)
         vy = my_r(ip1) / ro_r(ip1)
         vz = mz_r(ip1) / ro_r(ip1)
         b2 = bx_r(ip1) * bx_r(ip1) + by_r(ip1) * by_r(ip1) + bz_r(ip1) * bz_r(ip1)
         v2 = vx * vx + vy * vy + vz * vz
         pr = (en_r(ip1) - v2 * ro_r(ip1) / 2.0d0 - b2 / 2.0d0) * gmm1
         c2 = gm * pr
         s2 = c2 + b2
         ca2 = bx_r(ip1) * bx_r(ip1)
         cfx = sqrt ((s2 + sqrt (s2 * s2 - 4.0d0 * c2 * ca2)) / ro_r(ip1) / 2.0d0)
         max_speed = max (cfx + abs (vx), max_speed)

         fro_r = mx_r(ip1)
         fmx_r = vx * mx_r(ip1) + pr + b2 / 2.0d0 - bx_r(ip1) * bx_r(ip1)
         fmy_r = vx * my_r(ip1) - bx_r(ip1) * by_r(ip1)
         fmz_r = vx * mz_r(ip1) - bx_r(ip1) * bz_r(ip1)
         fbx_r = 0.0d0
         fby_r = vx * by_r(ip1) - bx_r(ip1) * vy
         fbz_r = vx * bz_r(ip1) - bx_r(ip1) * vz
         fen_r = (en_r(ip1) + pr + b2 / 2.0d0) * vx - (bx_r(ip1) * vx + by_r(ip1) * vy + bz_r(ip1) * vz) * bx_r(ip1)

         fro(i) = 0.5d0 * (fro_l + fro_r - max_speed * (ro_r(ip1) - ro_l(i)))
         fmx(i) = 0.5d0 * (fmx_l + fmx_r - max_speed * (mx_r(ip1) - mx_l(i)))
         fmy(i) = 0.5d0 * (fmy_l + fmy_r - max_speed * (my_r(ip1) - my_l(i)))
         fmz(i) = 0.5d0 * (fmz_l + fmz_r - max_speed * (mz_r(ip1) - mz_l(i)))
         fbx(i) = 0.5d0 * (fbx_l + fbx_r - max_speed * (bx_r(ip1) - bx_l(i)))
         fby(i) = 0.5d0 * (fby_l + fby_r - max_speed * (by_r(ip1) - by_l(i)))
         fbz(i) = 0.5d0 * (fbz_l + fbz_r - max_speed * (bz_r(ip1) - bz_l(i)))
         fen(i) = 0.5d0 * (fen_l + fen_r - max_speed * (en_r(ip1) - en_l(i)))
      enddo
   endif

   do i = 2, nx
      im1 = i - 1
      ro(i) = ro(i) - dtdx * (fro(i) - fro(im1))
      mx(i) = mx(i) - dtdx * (fmx(i) - fmx(im1))
      my(i) = my(i) - dtdx * (fmy(i) - fmy(im1))
      mz(i) = mz(i) - dtdx * (fmz(i) - fmz(im1))
      bx(i) = bx(i) - dtdx * (fbx(i) - fbx(im1))
      by(i) = by(i) - dtdx * (fby(i) - fby(im1))
      bz(i) = bz(i) - dtdx * (fbz(i) - fbz(im1))
      en(i) = en(i) - dtdx * (fen(i) - fen(im1)) 
   enddo

   flux_bnd(1, 1) = fro(ng)
   flux_bnd(2, 1) = fro(nx - ng)
   flux_bnd(1, 2) = fmx(ng)
   flux_bnd(2, 2) = fmx(nx - ng)
   flux_bnd(1, 3) = fmy(ng)
   flux_bnd(2, 3) = fmy(nx - ng)
   flux_bnd(1, 4) = fmz(ng)
   flux_bnd(2, 4) = fmz(nx - ng)
   flux_bnd(1, 5) = fbx(ng)
   flux_bnd(2, 5) = fbx(nx - ng)
   flux_bnd(1, 6) = fby(ng)
   flux_bnd(2, 6) = fby(nx - ng)
   flux_bnd(1, 7) = fbz(ng)
   flux_bnd(2, 7) = fbz(nx - ng)
   flux_bnd(1, 8) = fen(ng)
   flux_bnd(2, 8) = fen(nx - ng)

   return
end subroutine tvd_lf
